Comparison  of  a  Finite  Difference  and  a  Mixed  Finite  Element 
Formulation  of  the  Uniaxial  Perfectly  Matched  Layer 

V.  A.  Bokif  and  M.  W.  Buksasb 

Center  for  Research  in  Scientific  Computation3, 

North  Carolina  State  University 
Raleigh,  NC  27695-8205 
vabokil@ncsu.edu 

CCS-4,  Mail  Stop  D409b 
Los  Alamos  National  Laboratory 
Los  Alamos,  NM  87544 
mwbuksas@lanl.gov 

April  14,  2006 

Abstract:  We  consider  the  anisotropic  uniaxial  formulation  of  the  perfectly  matched  layer 
model  (UPML).  We  prove  the  decay  of  different  energies  for  the  UPML,  under  certain  as¬ 
sumptions,  to  demonstrate  the  well-posedness  of  this  formulation.  We  present  and  analyze 
a  mixed  finite  element  method  for  the  time  domain  discretization  of  the  UPML  to  simu¬ 
late  wave  propagation  in  unbounded  domains  in  two  dimensions.  On  rectangles  the  spatial 
discretization  uses  bilinear  finite  elements  for  the  electric  held  and  the  lowest  order  Raviart- 
Thomas  divergence  conforming  elements  for  the  magnetic  held.  We  use  a  centered  finite 
difference  method  for  the  time  discretization.  We  compare  the  finite  element  technique  pre¬ 
sented  to  the  finite  difference  time  domain  method  (FDTD)  via  a  stability,  dispersion,  phase 
error  and  numerical  reflection  coefficient  analysis.  We  derive  the  reflection  coefficient  for 
the  case  of  a  semi-infinite  layer  to  show  consistency  between  the  numerical  and  continuous 
models,  and  in  the  case  of  a  finite  PML  to  study  the  effects  of  terminating  the  absorbing 
layer.  Finally,  we  demonstrate  the  effectiveness  of  the  mixed  finite  element  scheme  by  nu¬ 
merical  examples  and  provide  comparisons  with  the  split  held  PML  discretized  by  the  FDTD 
method.  In  conclusion,  we  observe  that  the  mixed  finite  element  scheme  for  the  PML  model 
has  absorbing  properties  that  are  comparable  to  the  FDTD  method. 

Keywords:  Perfectly  matched  layers,  Mixed  finite  element  methods,  FDTD,  Maxwell’s 
equations. 
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1  Introduction 


The  effective  modeling  of  waves  on  unbounded  domains  by  numerical  methods,  such  as  the 
finite  difference  method  or  the  finite  element  method,  is  dependent  on  the  particular  ab¬ 
sorbing  boundary  condition  used  to  truncate  the  computational  domain.  In  1994,  J.  P. 
Berenger  created  the  perfectly  matched  layer  (PML)  technique  for  the  reflectionless  absorp¬ 
tion  of  electromagnetic  waves  in  the  time  domain  [6].  The  PML  is  an  absorbing  layer  that  is 
placed  around  the  computational  domain  of  interest  in  order  to  attenuate  outgoing  radiation. 
Berenger  showed  that  his  PML  model  allowed  perfect  transmission  of  electromagnetic  waves 
across  the  interface  of  the  computational  domain  regardless  of  the  frequency,  polarization 
or  angle  of  incidence  of  the  waves.  The  waves  are  then  attenuated  exponentially  with  depth 
in  the  absorbing  layers.  Since  its  original  inception  in  1994  PML’s  have  also  extended  their 
applicability  in  areas  other  than  computational  electromagnetics  such  as  acoustics,  elasticity, 
etc.  [2,  3,  17,  18,  19]. 

The  properties  of  the  continuous  PML  model  have  been  studied  extensively  and  are  well 
documented.  The  original  split  field  PML  proposed  by  Berenger  involved  a  nonphysical 
splitting  of  Maxwell’s  equations  resulting  in  non-Maxwellian  fields  and  a  weakly  hyper¬ 
bolic  system  [1],  A  complex  change  of  variables  approach  was  used  in  [11,  30]  to  derive  an 
equivalent  PML  model  that  did  not  require  a  splitting  of  Maxwell’s  equations.  In  [32],  the 
authors  observed  that  a  material  can  possess  reflectionless  properties  if  it  is  assumed  to  be 
anisotropic.  A  single  layer  in  this  technique  was  termed  uniaxial  and  the  PML  was  referred 
to  as  the  uniaxial  PML  (UPML).  In  this  method  modifications  to  Maxwell’s  equations  are 
also  not  required  and  one  obtains  a  strongly  hyperbolic  system.  In  [16,  26]  further  study  of 
the  anisotropic  PML  is  carried  out.  Unlike  Berenger ’s  split  field  PML,  which  is  a  nonphysical 
medium,  the  anisotropic  PML  can  be  a  physically  realizable  medium  [30].  Thus,  there  are 
several  reasons  for  using  the  anisotropic  PML  in  numerical  simulations.  In  [38]  the  authors 
show  that  the  anisotropic  PML  and  Berenger’s  split  field  PML  produce  the  same  tangential 
fields,  however,  the  normal  fields  are  different  as  the  two  methods  satisfy  different  divergence 
conditions. 

The  finite  depth  of  the  absorbing  layer  allows  the  transmitted  part  of  the  wave  to  return  to 
the  computational  domain.  In  addition,  the  discretization  of  Maxwell’s  equations  introduces 
errors  which  cause  the  PML  to  be  less  than  perfectly  matched.  Even  so,  it  has  been  found  that 
the  PML  medium  can  result  in  reflection  errors  as  minute  as  -80  dB  to  -100  dB  [6,  7,  11,  16]. 
There  are  a  number  of  publications  that  study  the  properties  of  the  finite  difference  time 
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domain  (FDTD)  method  (Yee  scheme  [41]),  applied  to  the  PML  model  (e.g.,  see  [33]).  There 
are  many  open  questions  related  to  the  numerical  approximations  to  the  PML  equations. 
For  example,  there  are  empirical  results  related  to  the  selection  of  the  PML  conductivities, 
but  no  rigorous  proofs  on  the  optimal  values  of  these  parameters  which  govern  the  numerical 
reflections  produced  due  to  discretization.  In  [14,  40,  28,  20]  the  authors  study  numerical 
errors  in  the  discretized  PML  and  consider  various  methods  for  the  selection  of  optimal 
values  of  the  PML  conductivities. 

There  are  significantly  less  publications  that  study  the  properties  of  the  finite  element 
method  for  the  approximation  of  the  PML  equations.  A  comparison  of  the  anisotropic  PML 
to  the  split  held  PML  of  Berenger  was  performed  in  [38] ,  in  which  the  authors  implement  the 
anisotropic  PML  into  an  edge  based  finite  element  method  for  a  second  order  formulation  of 
Maxwell’s  equations.  In  [39]  the  authors  use  the  lowest  order  as  well  as  first  order  tangential 
vector  finite  element  methods  for  the  discretization  of  the  electric  held.  They  compare  the 
performance  of  these  elements  with  the  FDTD  method  when  a  PML  is  used  to  terminate 
the  computational  domain.  They  show  that  the  lowest  order  elements  do  not  perform  as 
well  as  the  FDTD  method,  however,  the  hrst  order  elements  can  produce  more  accurate 
results  than  FDTD.  The  lowest  order  edge  element  is,  however,  a  most  commonly  used 
finite  element  approximation  in  the  literature  and  most  of  these  implementations  are  in  the 
frequency  domain  (for  e.g.,  see  [25,  34,  22]).  A  time  domain  mixed  finite  element  method 
has  been  used  in  [13]  along  with  mass  lumping  techniques  to  solve  scattering  problems  on 
domains  where  a  PML  method  based  on  the  Zhao-Cangellaris’s  model  is  used  to  terminate 
the  mesh  [42],  The  underlying  partial  differential  equations  in  the  Zhao-Cangellaris’s  PML 
model  are  second  order  in  time,  whereas  the  anisotropic  uniaxial  model  consists  of  a  system 
of  hrst  order  PDE’s. 

In  this  paper,  we  present  a  mixed  finite  element  method  (FEM)  for  the  discretization  of 
the  anisotropic  uniaxial  formulation  of  the  UPML,  by  Sacks  et  al,  [32]  in  the  time  domain 
to  simulate  wave  propagation  on  unbounded  regions.  We  divide  the  computational  domain 
into  rectangles.  On  each  rectangle,  we  use  continuous  piecewise  bilinear  finite  elements  to 
discretize  the  electric  held  and  the  Raviart-Thomas  elements  [31]  to  discretize  the  magnetic 
held.  The  degrees  of  freedom  are  staggered  in  space  as  in  the  FDTD  scheme.  We  use  a 
centered  difference  scheme  in  time  again  staggering  the  temporal  components  of  the  electric 
and  magnetic  helds.  We  study  the  numerical  properties  of  the  presented  method  and  the 
effectiveness  of  the  PML  technique  as  an  absorbing  boundary  condition  for  FEM.  We  provide 
comparisons  of  the  numerical  approximations  of  the  PML  model  by  the  mixed  FEM  with 
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those  by  the  FDTD  method.  The  presented  finite  element  method  can  be  viewed  as  an 
averaged  FDTD  scheme  on  regular  grids.  We  have  used  this  method  in  problems  of  scattering 
type  in  [9].  A  mixed  finite  element  method  was  used  for  the  discretization  of  a  similar  UPML 
formulation  for  the  wave  equation  written  as  a  system  of  first  order  PDE’s  in  [8]. 

An  advantage  of  the  FEM  is  that  it  can  model  arbitrary  complex  geometrical  structures 
effectively.  We  compare  the  UPML  discretized  by  the  finite  element  method  with  the  split 
held  PML  of  Berenger  discretized  by  the  FDTD  method.  These  comparisons  demonstrate 
that  the  PML  technique  is  an  effective  absorbing  boundary  condition  for  the  mixed  finite 
element  method  which  has  comparable  (with  FDTD)  absorbing  properties.  This  implies 
that  using  mass  lumping  techniques  [5]  one  can  use  the  PML  method  as  a  good  absorbing 
boundary  condition  with  FEM  for  problems  involving  complicated  geometrical  structures. 
The  mass  lumping  techniques  can  help  in  significantly  lowering  costs  associated  with  the 
implementation  of  the  finite  element  method  and,  thus,  make  it  a  more  general  alternative 
to  the  FDTD  method  for  discretizing  PML  models. 

In  Sections  2  and  3,  we  describe  the  UPML  model  and  its  implementation.  In  Section  4, 
we  derive  the  two-dimensional  (2D)  transverse  magnetic  (TM)  mode  of  the  UPML  model,  and 
we  describe  a  mixed  finite  element  formulation  for  the  LIPML.  We  state  some  energy  decay 
results  that  indicate  the  well-posedness  of  the  PML  model  in  Section  5.  Section  6  describes 
the  numerical  discretization  in  space  and  time.  We  perform  a  dispersion  and  stability  analysis 
in  Section  7.1,  a  phase  error  analysis  in  Section  7.2  and  a  reflection  coefficient  analysis  in 
Section  7.3.  In  both  these  analyses,  we  provide  comparison  of  the  properties  of  the  mixed 
finite  element  method  with  the  FDTD  method.  Finally,  we  present  numerical  examples  in 
Section  8  that  demonstrate  the  effectiveness  of  the  discrete  PML  model  for  FEM.  We  also 
compare  these  numerical  results  with  those  of  the  split  held  method  of  Berenger  discretized 
by  the  FDTD  method. 
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2  An  Anisotropic  Perfectly  Matched  Layer  Model 


We  begin  with  a  form  of  Maxwell’s  equations,  suitable  for  general  media,  which  permits  both 
electric  and  magnetic  currents  but  does  not  contain  unbalanced  electric  charges 

(9B 

-7 —  =  — VxE  —  Jm,  (Maxwell-Faraday’s  Law), 


<9D 

-jj-  =  VxH  —  Je,  (Maxwell- Ampere’s  Law), 

V  •  B  =  0,  (Gauss’s  Law  for  the  magnetic  field), 


(2.1) 


V  •  D  =  0,  (Gauss’s  Law  for  the  electric  field). 

Constitutive  relations  which  relate  the  electric  and  magnetic  fluxes  (D,B)  and  the  electric 
and  magnetic  currents  (Je,  Jm)  to  the  electric  and  magnetic  fields  (E.  H)  are  added  to 
these  equations  to  make  the  system  fully  determined  and  to  describe  the  response  of  a 
material  to  the  electromagnetic  fields.  In  free  space,  these  constitutive  relations  are  D  = 
eoE  and  B  =  /r0H,  and  Je  =  Jm  =  0,  where  eo,  and  /i0  are  the  permittivity,  and  the 
permeability  of  free  space,  respectively.  In  general,  there  are  different  possible  forms  for  these 
constitutive  relationships.  In  a  frequency  domain  formulation  of  Maxwell’s  equations,  these 
can  be  converted  to  linear  relationships  between  the  dependent  and  independent  quantities 
with  frequency  dependent  coefficient  parameters. 

We  will  derive  a  PML  model  in  the  frequency  domain  and  then  obtain  a  PML  model  in 
the  time  domain  by  taking  the  inverse  Fourier  transforms  of  the  frequency  domain  equations. 
To  this  end,  we  consider  the  time-harmonic  form  of  Maxwell’s  equations  (2.1)  (with  time 
dependence  elut)  given  by 

icuB  =  -VxE  -  JM, 
icjD  =  VxH  —  3  Ei 

(2.2) 

V  •  B  =  0, 

V  -  D  =  0, 

where  for  every  field  vector  V,  V  denotes  its  Fourier  transform  and  we  have  the  constitutive 
laws 

B  =  [/x]H,  D  =  [e]E.  Jm  =  [<Jm\ H,  J  e  =  [cte]E.  (2.3) 
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Here,  the  square  brackets  indicate  a  tensor  quantity.  Note  that  when  the  density  of  electric 
and  magnetic  charge  carriers  in  the  medium  is  uniform  throughout  space,  then  V  ■  Je  =  0 
and  V  •  JM  =  0.  We  define  new  tensors 


[/i]  =  \n\  + 


Wm] 


e  —  e  + 


We] 


(2.4) 


iuj  io ; 

Using  the  definitions  (2.4)  we  define  two  new  constitutive  laws  that  are  equivalent  to  (2.3), 
given  by 


Bnew  Dnew  p  E. 


(2.5) 


Using  (2.5)  in  (2.2)  Maxwell’s  equations  in  time-harmonic  form  become 


iooBnew  =  —  VxE, 


icuDnew  =  V  X  H, 

(2.6) 

V  •  Bnew  =  0, 

V  •  D„ew  =  0. 

The  split-held  PML  introduced  by  Berenger  [6]  is  a  hypothetical  medium  based  on  a 
mathematical  model.  In  [26]  Mittra  and  Pekel  showed  that  Berenger’s  PML  was  equivalent 
to  Maxwell’s  equations  with  a  diagonally  anisotropic  tensor  appearing  in  the  constitutive 
relations  for  D  and  B.  For  a  single  interface  the  anisotropic  medium  is  uniaxial  and  is 
composed  of  both  the  electric  permittivity  and  magnetic  permeability  tensors.  This  uniaxial 
formulation  performs  as  well  as  the  original  split-held  PML  while  avoiding  the  nonphysical 
held  splitting.  As  will  be  shown  below,  by  properly  defining  a  general  constitutive  tensor  [S], 
we  can  use  the  UPML  in  the  interior  working  volume  as  well  as  the  absorbing  layer.  This 
tensor  provides  a  lossless  isotropic  medium  in  the  primary  computation  zone,  and  individual 
UPML  absorbers  adjacent  to  the  outer  lattice  boundary  planes  for  mitigation  of  spurious 
wave  rehections.  The  helds  excited  within  the  UPML  are  also  plane  wave  in  nature  and 
satisfy  Maxwell’s  curl  equations. 

The  derivation  of  the  PML  properties  for  the  tensor  constitutive  laws  is  also  done  directly 
by  Sacks  et  al.,  in  [32]  and  by  Gedney  in  [16].  We  follow  the  derivation  by  Sacks  et  al,  here. 
We  begin  by  considering  planar  electromagnetic  waves  in  free  space  incident  upon  a  PML 
half  space.  Starting  with  the  impedance  matching  assumption,  i.e.,  the  impedance  of  the 
layer  must  match  that  of  free  space:  Vo  =  [e]  _ 1  [/i]  we  have 

—  =  —  =  [£]  =  diag{ai,  ai,  a3}.  (2.7) 

eo  ho 
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Hence,  the  constitutive  parameters  inside  the  PML  layer  are  [e]  =  eofS],  and  [p]  =  /io  [S'] , 
where  [S]  is  a  diagonal  tensor. 

We  consider  plane  wave  solutions  of  the  form  V(x,  t)  =  V(x)e1^_k'x\  for  all  held  vectors 
V,  to  the  time-harmonic  Maxwell’s  equations  with  the  diagonally  anisotropic  tensor  [S].  Here 
k  =  (kx,  ky,  kz)  is  the  wave  vector  of  the  planar  electromagnetic  wave  and  x  =  (x,  y,  z ).  In 
this  case,  the  dispersion  relation  for  waves  in  the  PML  are  found  to  be 

1,2  v.2  y2  ,2 

+  —  =  kl  =  cn2/i0e0  =  -T,  (2.8) 

$2^3  CL\&2t  CL\CL2  Cz 

where  c  is  the  speed  of  light  in  free  space. 

Without  loss  of  generality  we  consider  a  PML  layer  which  fills  the  positive  x  half-space 
and  plane  waves  with  wave  vectors  in  the  xy-  plane  (kz  =  0).  Let  9t  be  the  angle  of  incidence 
of  the  plane  wave  measured  from  the  normal  to  the  surface  x  —  0.  The  standard  phase  and 
magnitude  matching  arguments  at  the  interface  yield  a  generalization  of  Snell’s  law 


yh/ 1  a  3  sin  6t  =  sin  0t, 


(2.9) 


where  9t  is  the  angle  of  the  transmitted  plane  wave.  By  matching  the  magnitudes  of  the 
electric  and  magnetic  fields  at  the  interface,  x  —  0,  the  reflection  coefficients  for  the  transverse 
electric  (TE)  and  the  TM  modes  are  given  to  be 


(2.10) 


From  (2.10)  we  see  that  by  choosing  a3  =  a2  =  a  and  \/aia3  =  1,  the  interface  is  completely 
reflectionless  for  any  frequency  and  angle  of  incidence  and  polarization.  Using  (2.5)  and 
(2.7)  the  constitutive  laws  for  the  perfectly  matched  layer  are 


Bnew  ho[‘S']H,  Dnew  CofS^E, 


where  the  tensor  [S]  is  given  as 


[S\ 


a -1  0  0 
0  a  0 
0  0  a 


(2.11) 


(2.12) 


The  perfectly  matched  layer  is  therefore  characterized  by  the  single  complex  number  a  in 
the  definition  of  the  tensor  [S'].  Taking  it  to  be  the  constant  a  =  7  —  i[3,  and  substituting 
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into  the  dispersion  relation  (2.8),  we  can  obtain  the  expression 


E (x,y,  z )  =  E0e-kopcosetxe-iko^cosdtx+sin9ty)eiut,  (2.13) 

for  the  electric  held  inside  of  the  PML.  Hence,  we  can  see  that  7  determines  the  wavelength 
of  the  wave  in  the  PML,  and  for  /3  >  0  the  wave  is  attenuated  according  to  the  distance  of 
travel  in  the  x  direction. 

3  Implementation  of  the  Uniaxial  PML 

To  apply  the  perfectly  matched  layer  to  electromagnetic  computations  the  half  infinite  layer 
is  replaced  with  a  layer  of  finite  depth  and  backed  with  a  more  conventional  boundary 
condition,  such  as  a  perfect  electric  conductor  (PEC).  This  truncation  of  the  layer  will  lead 
to  reflections  generated  at  the  PEC  surface,  which  can  propagate  back  through  the  layer 
to  re-enter  the  computational  region.  In  this  case,  the  reflection  coefficient  R  is  a  function 
of  the  angle  of  incidence  6,  the  depth  of  the  PML  <5,  as  well  as  the  parameter  a  in  (2.12). 
Thus,  the  parameter  a  for  the  PML  is  chosen  in  order  for  the  attenuation  of  waves  in  the 
PML  to  be  sufficient  so  that  the  waves  striking  the  PEC  surface  are  negligible  in  magnitude. 
Perfectly  matched  layers  are  then  placed  near  each  edge  (face  in  3D)  of  the  computational 
domain  where  a  non-reflecting  condition  is  desired.  This  leads  to  overlapping  PML  regions 
in  the  corners  of  the  domain.  As  shown  in  [32],  the  correct  form  of  the  tensor  which  appears 
in  the  constitutive  laws  for  these  regions  is  the  product 

[S]  =  [S]«[S]»[S]„  (3.1) 

where  component  [S']Q  in  the  product  in  (3.1)  is  responsible  for  attenuation  in  the  a  direction, 
for  a  =  x,y,z ;  see  Figure  1.  All  three  of  the  component  tensors  in  (3.1)  are  diagonal  and 
have  the  forms 

s”1  0  0  sy  0  0  0  0 

[S\x=  0  sx  0  ;  [S\y=  0  s'1  0  ;  [S]z  =  0  sz  0  .  (3.2) 

0  0  sx  0  0  sy  0  0  S71 

In  the  above  sx,sy,sz  are  analogous  to  the  complex  valued  parameter  a  encountered  in 
Section  2,  in  the  analysis  of  a  single  PML  layer.  Here  sa  governs  the  attenuation  of  the 
electromagnetic  waves  in  the  a  direction  for  a  —  x,y,  z. 

When  designing  PML’s  for  implementation  it  is  important  to  choose  the  parameters  sa 
so  that  the  resulting  frequency  domain  equations  can  be  easily  converted  back  into  the  time 
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[s]  =  [«t  ,,  ,, 

■ _ _ _ _ _  r  oi  roi 

II  i 

crx  >  0 

llllllllllllllllll  1  M 

§ 

Gy 

>0  Gy  =  0  Gy 
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Figure  1:  PML  layers  surrounding  the  domain  of  interest.  The  PML  conductivities  ax  and  ay 
are  both  zero  in  the  working  volume.  The  PML  is  truncated  by  a  perfect  electric  conductor. 


domain.  The  simplest  of  these  which  we  employ  here  [16]  is 


sa  =  1  H - — ,  where  aa  >  0,  a  =  x,  y,  z. 

icue0 


(3.3) 


The  PML  interface  represents  a  discontinuity  in  the  conductivities  aa.  To  reduce  the  numer¬ 
ical  reflections  caused  by  these  discontinuous  conductivities  the  aa  are  chosen  to  be  functions 
of  the  variable  a  (for  e.g.,  gx  is  taken  to  be  a  function  of  x  in  the  [S]x  component  of  the 
PML  tensor).  Choosing  these  functions  so  that  cra  =  0,  i.e.,  sa  —  1  at  the  interface  makes 
the  PML  a  continuous  extension  of  the  medium  being  matched  and  reduces  numerical  reflec¬ 
tions  at  the  interface.  Increasing  the  value  of  cra  with  depth  in  the  layer  allows  for  greater 
overall  attenuation  while  keeping  down  the  numerical  reflections.  Gedney  [16]  suggests  a 
conductivity  profile 


<ja  (o) 


^max|^  ^0 


m 


8m 


a  =  x,  y,  z, 


(3.4) 


where  8  is  the  depth  of  the  layer,  a  =  cko  is  the  interface  between  the  PML  and  the  com¬ 
putational  domain,  and  m  is  the  order  of  the  polynomial  variation.  Gedney  remarks  that 
values  of  m  between  3  and  4  are  believed  to  be  optimal.  For  the  conductivity  profile  (3.4) 
the  PML  parameters  can  be  determined  for  given  values  of  m,  8,  and  the  desired  reflection 
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coefficient  at  normal  incidence  Rq  as 


(m  +  1)  ln(l/i?0) 

<7"'“  “  2 ZS  ’ 


(3.5) 


Z  being  the  characteristic  wave  impedance  of  the  PML.  Empirical  testing  suggests  [16]  that, 
for  a  broad  range  of  problems,  an  optimal  value  of  crmax  is  given  by 


m  +  1 

opt  lbQTtha^/e^'> 


(3.6) 


where  ha  is  the  space  increment  in  the  a  direction  and  er  is  the  relative  permittivity  of  the 
material  being  modeled.  In  the  case  of  free  space  the  relative  permittivity  is  er  —  1. 


4  A  Mixed  Finite  Element  Formulation  for  the  UPML 
in  Two  Dimensions 


From  the  time- harmonic  Maxwell’s  curl  equations  in  the  UPML  (2.6)  and  (2.11),  Ampere’s 
and  Faraday’s  laws  can  be  written  in  the  most  general  form  as 

(  kn/io[S']H  —  — V xE  5  (Maxwell-Faraday  s  Law), 


(4.1) 


[  icaeofS'jE  =  VxH  ;  (Maxwell- Ampere’s  Law). 


In  (4.1),  [S']  is  the  diagonal  tensor  defined  via  (3.1)-(3.5).  In  the  presence  of  this  diagonal 
tensor,  a  plane  wave  is  purely  transmitted  into  the  uniaxial  medium.  The  tensor  [S]  is  no 
longer  uniaxial  by  strict  definition  but  rather  is  anisotropic.  However,  the  anisotropic  PML 
is  still  referenced  as  uniaxial  since  it  is  uniaxial  in  the  non  overlapping  PML  regions. 

To  obtain  the  2D  model  of  the  UPML  we  assume  no  variation  in  the  z  direction  (i.e., 
dz  =  0).  In  the  2D  TM  mode  the  electromagnetic  field  has  three  components  EZ,HX,  and 
Hy.  Let  dq  =  dq  denote  the  derivative  w.r.t  g,  for  q  =  x,y,z,t.  In  this  case,  we  have  az  =  0 
and  sz  =  1  in  the  UPML,  and  the  time-harmonic  Maxwell’s  equations  (4.1)  in  the  uniaxial 
medium  can  be  written  in  scalar  form  as 

ica/io"  Hx  9yEz , 

sx 

kjfio-Hy  =  -dxEzy  (4.2) 

Sy 

,  i^^O  SXSyEZ  DXHy  dyHx. 

To  avoid  a  computationally  intensive  implementation  we  do  not  insert  the  expressions  for 
sx,  sy  and  sz  obtained  via  (3.3)  into  (4.2)  and  transform  to  the  time  domain.  Instead  we  define 
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suitable  constitutive  relationships  that  facilitate  the  decoupling  of  the  frequency  dependent 
terms  [33].  To  this  end,  we  introduce  the  fields 


Bx  Hxi 

By  =  Po  Sy'Hy,  (4.3) 

k  Dz  P'0  SyEz. 

Substituting  the  definitions  (4.3)  in  (4.2),  using  the  defining  relations  for  sx  and  sy  from 
(3.3),  and  then  transforming  into  the  time  domain  by  using  the  inverse  Fourier  transform 
yields  an  equivalent  system  of  time-domain  differential  equations  given  as 


1) 

dtBx 

=  -^Bx-dyEz, 
eo 

ii) 

dtHx 

—  dfBx  + 

Po  e0po 

iii) 

dtBy 

= - By  +  dxEz , 

eo 

IV) 

dtHy 

—  dtBy  +  By , 

P  0  eoPo 

v) 

dtDz 

=  Dz  +  dx  Hy  - 

eo 

Vi) 

d \EZ 

= - cyEz  H - d tD 

eo  e0 

We  can  rewrite  (4.4)  in  vector  form  as 


dtB  = 

--e2b 

eo 

— 

curl  E, 

dtU  = 

-dtB 

Po 

+ 

1  SiB 
eoPo 

d \D  = 

e0 

+ 

curl  H, 

dtE  = 

--OyE 

e0 

+ 

-dtD. 

eo 

In  the  above  H  =  (Hx,  Hy ),  B  =  (Bx,  By ),  E  =  Ez  and  D  =  D.  In  (4.5) 

Si  =  (  0  )  i  Sid"’  ’)■ 

U  U  7# 


s2=("»  0 

V  0  Gr 


Let  D  denote  the  computational  domain  in  R2.  We  denote  the  domain  D  including  the 
surrounding  finite  PML  layers  by  Q.  In  (4.5),  the  operators  denoted  by  curl,  and  curl  are 
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linear  differential  operators  which  are  defined  as 

curl  0  =  (dy<j>,  -dx<j>),  curl  v  =  dxvy  -  dyvx ,  V  v  =  (vx,vy)  G  X>'(fi)2,  V  0  G  V\Vt).  (4.7) 

Here  D'(TL)  is  the  space  of  distributions  on  fh  The  operator  curl  appears  as  the  (formal) 
transpose  of  the  operator  curl  [15],  i.e., 

(curl  v,(j>)  =  (v,curl0),  Vv  G  V\Q)2,  0  G  (4.8) 

Thus,  the  PML  model  consists  in  solving  system  (4.4)  (or  (4.5))  for  the  six  variables  Bx, 
By,  Hx,  Hy,  Dz,  Ez  in  Q,  with  PEC  conditions  on  (El  to  terminate  the  PML;  namely, 
n  x  E  =  0,  on  dfl,  where  n  is  the  outward  unit  normal  to  <90.  In  the  case  of  the  2D 
TM  mode  the  PEC  condition  translates  to  E  —  Ez  —  0,  on  dVt.  We  also  have  the  initial 
conditions 


E(x,0)  =  Eq,  D(x,0)  =  E0,  H(a;,  0)  =  H0,  B(a;,  0)  =  H0,  for  iGfl.  (4.9) 

We  consider  the  following  variational  formulation  of  system  (4.5)  which  is  suitable  for  dis¬ 
cretization  by  finite  elements. 


Find  (£(-,t),D(-,i),H(-,t),B(-,t))  G  H^(Q)  x  H^(Q)  x  [L2(fl)]2  x  [L2(fl)]2  such  that  for 
all  ^  G  [L2(D)]2,  for  all  0  G 


(  d 


—  /  B  ■  \I>  dx  = - /  IbB  •  \1>  dx 


dt 

v/H^dx  =  —  #-  /  B4dx  + 


eo  Jn 
1  d 


curl  E  ■  tL  dx, 


y(ig  dt 


in 


Jn 


ExB  •  dx, 


dt 

—  [  D  ■  (j)  dx  = - [  axD  ■  (f)  dx  +  [  curl  0  •  H  dx, 

dt 


d 


eo  Jn 
1 


1  d 


(  dt 


—  E  ■  6  dx  = - /  avE  ■  6  dx  -( - —  /  D  ■  6  dx. 


eo  do 


eg  dt 


(4.10) 


We  assume  that  the  fields  (E,  D.  H,  B)  are  sufficiently  differentiable  in  time.  We  note  that, 
for  E  G  L2(fi),  curl  E  =  (dyE,  —dxE)  G  [L2(f2)]2  implies  that  both  the  partial  derivatives  of 
E  must  be  in  L2(f2).  Hence  we  must  have  E  G  Hl{ fl). 


5  Energy  Estimates  for  the  UPML 

Maxwell’s  equations  form  a  symmetric  hyperbolic  system.  The  solution  to  such  a  system  is 
strongly  well-posed  [21].  It  is  natural  then  to  consider  the  well-posedness  of  the  different  PML 
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models.  Abarbanel  and  Gottlieb  [1]  showed  that  Berenger’s  split-field  PML  is  a  weakly  well- 
posed  system;  thus  instabilities  could  appear  in  numerical  implementations  of  this  model. 
In  [29]  the  authors  demonstrate  that  the  Zhao-Cangellaris’s  model  for  the  PML  is  strongly 
well-posed.  Becache  and  Joly  [4]  show  this  wcll-posedness  explicitly  by  presenting  energy 
decay  results  for  the  2D  TE  mode  of  this  model. 

We  derive  energy  decay  results  for  the  2D  TM  mode  of  the  UPML  in  two  cases;  a  a 
positive  constant  and  cr  e  L°°(fl).  We  have  derived  estimates  for  the  UPML  model  under 
the  same  conditions  as  done  in  [4]  for  the  Zhao-Cangellaris  model.  The  Zhao-Cangellaris’s 
model  established  the  equivalence  between  the  Chew-Weedon’s  PML  model  [11]  based  on 
coordinate  stretching  and  the  anisotropic  model  by  Sacks  et  al.  As  a  consequence,  the  energy 
decay  results  for  the  Zhao-Cangellaris’s  PML  and  the  UPML  appear  to  be  similar.  This  is 
to  be  understood  in  the  sense  that  the  definitions  of  the  energies  involved  are  identical  in 
the  second,  third,  and  fourth  estimate  and  almost  identical  in  the  first  estimate. 

To  simplify  the  analysis,  we  assume  that  €q  =  /j>q  =  1  in  the  rest  of  this  section.  We  also 
assume  that  we  start  with  zero  initial  conditions  in  (4.9).  Let  (•,  •)  denote  the  L2(D)  inner 
product  and  j|  •  1 1  L2(f2)  the  corresponding  norm. 

Energy  Estimate  1  Let  us  assume  that  the  computational  domain  D  has  a  PML  in  the 
region  x  >  0.  In  this  case,  ax  =  a  and  cry  =  0.  For  a  positive  constant  value  of  a,  the  energy 
£f  of  the  PML  system  defined  as 

£i  =  2  (H^IIl^q)  +  \\^tHy\\i?{n)  +  WaBy\\iP(ii)  +  II (^  +  <T) -^IIl2^))  >  (5-1) 

is  a  decreasing  function  of  time.  It  satisfies  the  identity 

fs^-2<r\\dtHt\\lHn)<0.  (5.2) 

Proof.  Consider  (4.4),  with  ax  =  a  and  cry  =  0.  Applying  the  operator  ( dt  +  cr)  to  (4.4,  i), 
the  operator  dt  to  (4.4,  ii)  and  combining  the  two  results,  we  get 

d?Hx  =  -dy(dt  +  a)E.  (5.3) 

We  note  that,  since  a  is  a  constant,  the  operator  ( dt  +  cr)  commutes  with  the  operators  dy 
and  dx.  Taking  the  inner  product  of  both  sides  of  (5.3)  with  d tHx  we  have 

\f  ll9.ff.llLs!)  =  -  (9»  (9.  +  °)  E,  S.ffJ  ■  (5.4) 
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Next,  apply  the  operator  ( dt  +  a)  to  (4.4,  iii)  to  get 

(d t  +  cr)  By  =  (dt  +  c)  dxE.  (5-5) 

Taking  inner  products  on  both  sides  of  (5.5)  with  dtBy  we  get 

(%By,  dtBy)  +  2^  (dtBy,  dtBy)  +  ^  (By,  dtBy)  =  ((dt  +  G)  d; XE . ,  dtBy )  -  (5.6) 

From  (4.4,  iv),  using  dtBy  =  dtHy  in  (5.6),  we  have 

2dt  +  2(7  (dtHy,dtHy)  +  \\<jBy\ £2^)  =  (dx  (dt  +  g)  E,  dtHy ) .  (5.7) 

Finally,  to  eliminate  Z),  we  apply  the  operator  dt  to  (4.4,  v)  and  the  operator  ( dt  +  g)  to 


(4.4,  vi),  combining  the  results  to  get 

(dt  +  a)  dtE  =  dtdxHy  -  dtdyHx.  (5.8) 

Taking  the  inner  products  of  both  sides  of  (5.8)  with  (dt  +  g)  E,  we  have 

2dt  +  ^Hl2(Q)  =  ( dtdxHy ,  (dt  +  a)  E)  —  ( dtdyHx ,  (dt  +  a)  E) .  (5.9) 

Integrating  by  parts  in  the  right  hand  side  of  (5.9)  we  have 

2dt  +  ^Hl2^)  =  _  (dtHy,  dx  (dt  +  a)  E)  +  (dtHx,  dy  (dt  +  g)  E) .  (5.10) 

Adding  (5.4),  (5.7)  and  (5.10)  and  using  definition  (5.1)  of  £f  we  have 

pf  +  2a\\dtH„\\lm  =  0.  (5.11) 


As  cr  >  0  we  obtain  (5.2)  from  (5.11).  This  implies  that  £f  is  a  decreasing  function  of  time. 

■ 

Energy  Estimate  2  Let  the  computational  domain  D  be  surrounded  by  PML ’s  on  all  sides 
( see  Figure  1).  Consider  a  corner  region  of  two  overlapping  layers  in  which  both  crx  and 
ay  are  positive  and  constants.  The  solution  of  the  UP  ML  TM  mode  satisfies  the  energy 
inequality 

S^ft)  <  £f(s),  for  all  t  >  s,  (5-12) 

where  * s  the  second  order  energy  defined  as 

=  l  (li3,2H||^(n)  +  ||E291H||J2(n)  +  ||(5,  +  (5,  +  a,)  £||2,(n))  .  (5.13) 
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Proof.  The  UPML  model  for  this  case  is 


(i)  dt  B  +  £2B  = 

—  curl  E, 

(ii)  dt  B  +  £iB  = 

dt  H, 

iii)  (dt  +  crx)  D  = 

curl  H. 

iv)  (dt  T  (Jy)  E  = 

dtD. 

(5.14) 


We  eliminate  B  from  (5.14,  i)  and  (5.14,  ii)  by  applying  the  operator  (dt  +  £2)  to  (5.14,  ii) 
and  the  operator  (dt  +  Id)  to  (5.14,  i)  and  combine  the  results  to  get 

dt(d*  +  £2)H  =  *-(dt  +  ?:1)wlE.  (5.15) 


Applying  the  operator  (dt  +  £2)  to  both  sides  in  the  equation  above,  we  get 

dt  ( dt  +  £2)2  H  =  -  (dt  +  £2)  (dt  +  £0  cVl  E.  (5.16) 


Taking  the  inner  product  of  both  sides  with  (9/H  we  have 

(dt  (dt  +  £2)2  H.  d; H)  =  -  ((dt  +  £2)  (dt  +  £1)  cVl  E,  d2t  h)  .  (5.17) 


This  implies 


=  -  ((0,  +  £2)  (Si  +  £1)  cud  E,  S2h) 


(5.18) 


Next,  to  eliminate  D  we  apply  the  operator  dt  to  (5.14,  iii),  and  the  operator  (dt  +  <JX ) 
to  (5.14,  iv),  and  combine  the  results  to  get 


(dt  +  crx)  (dt  +  a.y)  E  =  dt  curl  H.  (5.19) 

Applying  dt  to  both  sides  of  the  above  equation  and  taking  inner  products  with 
(dt  +  cfx)  (dt  +  (Jy)  E,  we  get 

(dt  (dt.  +  crx )  (dt  +  cry )  E,  (dt  +  ux)  (dt  +  <ry )  E)  =  (• d 2  curl  H.  (dt  +  ox )  (dt  +  cry )  A)  .  (5.20) 

This  implies 

II  (dt  +  °x)  (dt  +  <Jy)  E\\l2[n)  =  (d2t  curl  H,  (dt  +  <jx)  (dt  +  ay)  E)  .  (5.21) 
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Integrating  the  right  hand  side  by  parts  and  making  use  of  the  fact  that  the  operators 
(dt  +  crx )  and  ( dt  +  cry)  commute,  we  have 

1  d 

25ll(a«  +  <Tj(9,  +  a„)£fe(n|=  (9«2H,(9t  +  E2)(9,  +  E1)53E).  (5,22) 

Adding  (5.18)  and  (5.22),  using  the  definition  of  the  energy  £!f ,  and  as  crx  >  0  and  ay  >  0, 
we  have 

=  -2E2  ||3,2H|S^(n)  <  0,  (5.23) 

which  implies  that  Sfj''  is  a  decreasing  function  of  time.  Thus,  (5.12)  is  proved.  ■ 


Energy  Estimate  3  Again,  assume  a  PML  in  the  region  x  >  0.  For  a  G  L°°(Q),  the  zero 
order  energy  £f  of  the  UPML  defined  by 

M*)  =  2  (JI-^IIl2^)  +  ll-^2/llL2(n)  PIIl2(U)  +  W^x  II L2(r2)^)  j  (^-24) 

satisfies  the  energy  estimate 

SS(t)<£S(0)  +  2\\a\\oo  [t£S(s)ds.  (5.25) 

Jo 

Proof.  Consider  (4.4),  with  ox  =  a  and  ay  =  0.  Taking  the  inner  product  of  both  sides  of 
(4.4,  i)  with  HXl  both  sides  of  (4.4,  iii)  with  Hy ,  both  sides  of  (4.4,  v)  with  E,  adding  the 
three  resulting  equations  and  integrating  the  right  hand  side  by  parts  (IBP)  we  get 


(dtBx,  Hx)  +  ({dt  +  a)By,  Hy)  +  ((dt  +  a)D ,  E)  = 

-  (dyE,  Hx)  +  (dxE,  Hy)  +  (dxHy  -  8yHx,  E)  =  0  (by  IBP). 

Assuming  that  we  start  from  zero  initial  conditions,  from  (4.4,  ii,iv,  vi),  we  have 

(i)  Hx(',t)  =  Bx(-,t)  +  aBx(-,t), 

(ii)  Hy(-,t)  =  By(-,t), 

(iii)  E(-,t)  =  D(-,t), 


where,  in  the  above 


From  (5.27,  ii),  we  have 


Bx(-,t)  =  /  Bx(-,s )  ds. 


((dt  +  <j)By,  Hy)  =  (, dtHy ,  Hy)  +  (<7  Hy  ,  H y) 


1  d 

2  dt 


I H 


yiiL2(n) 


+  ( CHy ,  Hy)  . 


(5.26) 


(5.27) 


(5.28) 


(5.29) 
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Next,  from  (5.27,  iii),  we  have 


m  +  <7)0,  E)  =  (StE,  E)  +  (<70,  E)  =  ~  ||0|&(n)  +  (<70, 0) .  (5.30) 


Finally,  from  (4.4,  ii),  we  have 

(dtBx,  Hx)  =  (dtHx  -  aBx,  Hx)  =  (dtHx,  Hx)  -  (aBx,  Hx) . 

Using  (5.27,  i),  we  obtain 

(d,B„Hx)  =  \jt  ||0,|k(n)  -  (<7(0,  -oBx),Hx) 

=  \jt  l|0,||^(n|  -  (<70„0,)  +  (<72B„  ft)  . 


(5.31) 


(5.32) 


The  last  term  in  (5.32)  can  be  rewritten  using  (5.27,  i),  (4.4,  ii)  and  rearranging  terms  as 

a2Bx,Hx )  =  (, cr(Hx  -  BX),HX) 


=  ( a(Hx  -  BX),HX  -  Bx)  +  (dt(Hx  -  Bx),  Bx  -  Bx) . 


(5.33) 


Using  (5.33)  in  (5.32)  we  get 

(d,B„  H,)  =  ~  (||ftlk(n)  +  lift  -  s,||y(n)) 

-  (c tHx ,  Hx)  +  (a(Hx  -  BX),HX  -  Bx) . 

Substituting  (5.29),  (5.30)  and  (5.34)  in  (5.26)  we  have 

2~dt  H^2/llL2(n)  +  II ^IIl2(£7)  +  W^x  —  -®z|Il2(Q)) 

=  ( crHx ,  Hx)  -  ( oHy ,  By)  -  (crE,  E )  -  ( a(Hx  -  Bx),  Hx  -  Bx). 


(5.34) 


(5.35) 


The  right  hand  side  of  (5.35)  can  be  bounded  as 

(< oHx ,  Hx)  -  ( oHy ,  Hy)  -  ( aE ,  E)  -  ( a(Bx  -  Bx),  Bx  -  Bx) 

—  2  II  O'!!  oo  {  2  (|I^IIl2(Q)  +  II-^J/IIl2(0)  +  I|£|Il2(Q)  +  \\H-x  —  II  L2(r2)^)  )  =  2  Halloo  £o  ) 

where,  Sq  is  defined  in  (5.24).  Thus,  from  (5.24),  (5.35)  and  (5.36)  we  have 

Jt£ S(t)  <  2  IkIL  £„*(*)  =>  £S(t)  <  £„*(0)  +  2  I  ML  /'  £?(>)  ds.  (5.37) 
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Energy  Estimate  4  Let  the  computational  domain  D  be  surrounded  by  PML ’s  on  all  sides 
( see  Figure  1).  Consider  a  corner  region  of  two  overlapping  layers  in  which  both  ax  and  ay 
are  in  L°°(Q).  If  the  product  axay  remains  positive  everywhere  in  the  domain  of  interest 
then  the  zero  order  energy  of  the  UPML  given  by 

£o(t)  =  2  (j|H  —  B||L2(n)  +  1 1  1 1  l2  (f2)  +  l|H||L2(n)  +  (axcryE,  ,  (5.38) 


where  Eft )  =  /  Ef,s)ds,  satisfies  the  energy  estimate 


£o\t)  <  £o  (0)  +  3  (llcrariloo  +  INIJ  f  So  ( s)ds .  (5.39) 

Jo 

Proof.  The  UPML  equations  in  this  case  are  given  by  (5.14).  From  (5.14,  i,  iii),  we  get  after 
integrating  by  parts, 

((dt  +  E2)B,  H)  +  ((dt  +  ax)D,  E)  =  —  ^curl  E,  h)  +  (curl  H,  E)  =  0.  (5.40) 

Consider  the  term 

((dt  +  crx)D,  E)  =  (dtD,  E)  +  (axD,  E) .  (5.41) 

From  (5.14,  iv),  assuming  zero  initial  conditions  we  have 

D(-,t)  =  E(-,t)  +  <jyE(-,t).  (5-42) 


Thus,  from  (5.41)  and  (5.42)  we  have 

(( dt  +  <JX)D ,  E)  =  (dtE,  E)  +  (< txE ,  E)  +  (. cyE ,  E)  +  ^ crxayE ,  E^j 

=  2  dt  +  2dt  {axCryE,E)  +  ((^  +  ay)E'E)  • 

From,  (5.14,  ii),  we  have 

H(.,t)  =  B(-,t)  +  E1B(.,t), 


where 


B(-,t)=  /  B(-,  s)  ds 


From,  (5.44)  and  (5.45)  we  get 

((dt  +  £2)B,  H)  =  (3*H,  H)  +  ((ay  -  ax)Bx,  Hx)  +  ((ax  -  ay)By,  Hy 


1  d 

2  dt 


I H I 


L2(Q) 


+  (E2B,H)-(E1B,H) 


(5.43) 


(5.44) 

(5.45) 


(5.46) 
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Again,  from  (5.44)  and  (5.45)  we  get 


(SiB.  H)  =  (, axBx ,  Hx)  +  (. ayBy ,  Hy) 

^  (cTy(Hy  CTyBy^j  Hyj 

=  i>,  HJ:.  Hx)  +  (c lyHy,  Hy)  -  (  .t;  /;,  ,  //,  )  -  (crjB„,  //.,  ) 

=  (£,H,H)  -  (a2xB„Hx')  -  (^B„,B„)  . 

Using  (5.47)  in  (5.46)  we  get 

((9,  +  E2)B,  H)  =  ~  I|H||l2(i1,  +  (£2B,  H)  -  (E,H,  H) 

+  +  (a;By,Hy) 

=  \jt  l|Hfe(n)  +  (E2B,  H)  -  (E,H,  H)  +  (E?B.  H)  . 
We  can  simplify  the  last  term  in  (5.48)  as  follows 


(5.47) 


(5.48) 


(e?B,  h)  =  (Si(H  -  B),  H  -  B)  +  (SiH,  B)  +  ($(H  -  B),  H  -  B)  -  (E^,  H) 
=  \jt  ||H  -  B||^(n)  +  (E,(H  -  B),  H  -  B) . 

Also 

(S2B,  H)  =  (E2(B  -  H).  H)  +  (E2H,  H) . 

From  (5.48),  (5.49),  and  (5.50)  we  have 

((9,  +  E2)B,  H)  =  ~  (||H||^(0)  +  ||H  -  B||j>(n))  +  (E2(B  -  H),  H) 

+  ((E2  -  EOH,  H)  +  (E,(H  -  B),  H  -  B) 


(5.49) 


(5.50) 


(5.51) 


From  (5.43)  and  (5.51)  we  have 


m  +  <tx)D,  E)  +  ({dt  +  E2)B,  H)  =  (EUH  -  B),  H  -  B) 

+  2dt  (11^11^2(0)  +  HH!lL2(n)  +  llH  _  BIIl2(q)  +  fayE,E)) 

+  ((ax  +  ay)E,  E )  +  (E2(B  -  H).  H)  +  ((S2  -  E^H.  H)  =  0.  (5.52) 
Using  the  definition  of  the  energy  Sq  in  (5.38)  we  have 

jfS  =  ((Si  -  E2)H,  H)  -  (E,(H  -  B),  H  -  B)  -  ((<7*  +  a,)E,  E ) 

+  (E2(H-B),H). 
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We  can  bound  the  terms  on  the  right  hand  hand  of  (5.53)  as  follows.  We  have 


(E2(H-B),H)  +  ((E1-E2)H,H) 

<  (Moo  +  KlloJ  2  (|IHHl2(o)  +  llH  ~  bIIl2(0))  +  l|H||L2(n) 

Substituting  (5.54)  in  (5.53)  we  obtain 

jt£S  <  {IKIL  +  IWU  [iWlhn)  +  \  l|H||^(n|  +  \  IIH  -  B||=a(0) 

<  H\W*\L  +  \Wv\L}tf- 

Integrating,  we  thus  have  the  energy  estimate 

£o\t)  <  £o  (0)  +  3  (II^H^  +  ||aj  J  [  £$ (. s)ds , 

Jo 

which  is  (5.39).  ■ 

6  The  Discrete  Mixed  Finite  Element  Scheme 

6.1  Space  Discretization 

Let  O  be  a  union  of  rectangles  defining  a  regular  mesh  (7^)  with  square  elements  (K)  of 
edge  h  >  0  as  in  Figure  2.  We  consider  the  following  approximation  space  for  H  and  B: 

Vh  =  {*h  e  [L2m2 1  WK  e  Th,  Vh\K  e  RT[0]},  (6.1) 

where,  RT[0 ]  =  Pw  x  P0 1;  is  the  lowest  order  Raviart-Thomas  space  [31]  and  for  k i ,  k2  G 
NU{0}, 

Pk1k2  =  {p{xi,X2)\p(xUX2)  =  ^  Y  aVXlXl}- 

0<i<k!  0<j<k2 

The  basis  functions  for  Hx  have  unity  value  along  an  ey  and  are  zero  over  all  other  edges 
(see  Figure  2).  Similarly,  the  basis  functions  for  Hy  have  unity  value  along  an  ex  edge  and 
are  zero  over  all  other  edges  (see  Figure  2). 

The  approximation  space  for  E  and  D  is  chosen  to  be 

Uh  =  {ct>h  G  Hl(V) |  \/K  G  %Mk  G  Qi},  (6.2) 

where  the  space  Q\  —  P\\.  The  basis  functions  for  E  have  unity  value  at  one  node  and  are 
zero  at  all  other  nodes.  Figure  2  shows  the  locations  for  the  degrees  of  freedom  for  for  the 
electric  and  magnetic  fields. 


(5.54) 


(5.55) 


(5.56) 
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Edge  2 


E 

d 

g 

e 

3 


Edge  1 


Figure  2:  A  sample  domain  element  K .  The  degrees  of  freedom  for  the  electric  and 
magnetic  fields  are  staggered  in  space.  The  degrees  of  freedom  for  the  electric  field 
E  are  at  the  nodes  of  the  square.  The  degrees  of  freedom  for  Hx  and  Hy  are  at  the 
midpoints  of  edges  parallel  to  the  x-axis  and  y-axis,  respectively. 


Based  on  the  approximation  spaces  described  above  the  space  discrete  scheme  reads: 
Find  (Eh(-,t),Dh(-,t),Hh(-,t),Bh(-,t))  EUhxUhxVhx  Vh  such  that  for  all  \Eqt  e  Vh, 
for  all  4>h  E  Uh, 


C  —  f  B/j  •  'i'hdx.  = 
dt  Jo 
d  r 

—  Hh-  T'fedx  = 
dt  Jn 


h 


J  I  D 

dt  Jo 

i  UE* 


dx  = 
dx  = 


- /  EiBh  •  T'fedx  - 

eo  Jo 

1  d  /‘  T  , 

—  —  /  Bh  ■'S'hdx.  + 

do  vt  Jo 

[  VxDh  ■  (fih  dx  + 
eo  Jo 

1 


curl  Eh  ■  fl/ftdx, 


1 


E2Bh  ■  M>ftdx, 

e0d0  Jo 

/  curl  J)h  ■  Hh  dx, 
Jo 
1  d 


(6.3) 


- /  VyEh  '  < hdx  +  —  -r,  \  Dh  ■  4>h  dx 


eo  Jo 


eo  dt 


6.2  Time  Discretization 


For  the  time  discretization  we  use  a  centered  second  order  accurate  finite  difference  scheme. 
For  k  E  Z  let, 


yk+l/2  _  yk-1/2  yk+1/2  ,  yk- 1/2 

D^‘  = - A( - ■  ^  1 - ■  <6-4) 
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Using  the  definitions  (6.4)  we  can  describe  the  fully  discrete  scheme  in  space  and  time  as: 

Find  (E%+\D%+\  Hnh+\Bnh+^)  G  lAh  x  lAh  x  Vh  x  Vh  such  that  for  all  Vh  e  Vh,  for  all 
(, t>h  £  Uh-, 


'  (i) 

(DAtBnh,Vh)  = 

b) 

1  , 

— 

(cMK  , 

1  ,  _ 

(ii) 

(DAtUl^h)  = 

~(DAtBnh,*h) 

p  o 

+ 

(S1 

to/io 

(iii) 

(- DAtDnh+1\4>h )  = 

+ 

(curl  (j>h,  H"+"), 

(iv) 

(daX+Km  = 

-  —  (ayEhn+2  ,(j)h) 
e0 

+ 

~(DAtDnh+\  <j>h) 
Co 

7  Analysis  of  the  Discrete  PML  Model 

The  numerical  approximation  of  time-dependent  wave  problems  introduces  errors  which 
involve  dissipation ,  dispersion,  and  anisotropy.  The  attenuation  of  the  amplitude  of  the  plane 
wave  is  referred  to  as  dissipation.  The  mixed  finite  element  scheme  presented  in  Section  5.7 
is  a  non-dissipative  scheme  for  Maxwell’s  equations  in  free  space  [24].  However,  the  PML 
model,  in  which  lower  order  terms  are  added,  is  a  dissipative  model.  As  seen  in  (2.13)  waves 
are  attenuated  in  the  PML  according  to  the  distance  of  travel  in  a  given  direction. 

The  numerical  model  produces  waves  that  propagate  at  incorrect  speeds.  The  dependence 
of  the  velocity  of  propagation  of  the  numerical  sinusoidal  waves  on  frequency  is  termed  as 
dispersion,  while  the  dependence  of  the  velocity  upon  direction  is  referred  to  as  numeri¬ 
cal  anisotropy.  These  errors  are  cumulative,  which  implies  that,  as  waves  propagate  over 
long  distances  the  numerical  solution  becomes  corrupted  and  may  completely  deviate  from 
the  correct  solution.  A  dispersive  equation  is  one  that  admits  plane  wave  solutions  of  the 
form  el(w<~k  for  which  the  speed  of  propagation  is  dependent  on  the  frequency.  For  time- 
harmonic  waves,  numerical  dispersion  results  in  the  creation  of  a  phase  error  in  the  solution. 
This  is  due  to  the  incorrect  modeling  of  the  sinusoidal  behavior  of  the  propagating  wave,  as 
the  piecewise  polynomial  approximation  of  a  finite  element  method  does  not  exactly  match 
a  sine  or  cosine  function  [23]. 

In  this  section  we  study  the  properties  of  the  discrete  PML  model  in  terms  of  a  numerical 
dispersion  analysis  which  is  performed  on  regular  square  elements.  It  describes  the  propaga¬ 
tion  of  waves  in  the  numerical  method  away  from  the  boundaries.  In  addition,  we  also  obtain 
information  about  the  expected  accuracy  of  the  method.  We  also  compare  the  stability  and 
dispersion  properties  of  the  finite  element  method  with  those  of  the  finite  difference  method 
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A  study  of  the  dispersion  analysis  of  different  mixed  methods  for  the  time  domain  Maxwell’s 
equations  is  done  in  [27]. 

Recall  that  k  denotes  the  wave  vector  for  the  continuous  case.  For  simplicity,  we  again 
assume  that  e0  =  Po  =  1,  hence  c—1.  We  look  for  solutions  to  the  continuous  system  (4.10), 
of  the  form 

V(x,y,t)  =  V0e™-k~  (7.1) 

where  V  is  any  one  of  E,  D,  Hx,  Hy,  Bx,  By.  Substituting  (7.1)  in  (4.10)  shows  that  oo  and 
k  are  related  by  the  dispersion  relation 


oo 


2 


(7.2) 


Inside  the  computational  domain,  where  sx  =  sy  —  1,  the  dispersion  relation  is  given  by 


u?  =  k: 


+  K  = 


oo  =  |k| 


(7.3) 


Other  solutions  are  oo  —  0,  or  oo  —  -*-  |k|.  There  are  two  types  of  velocities  that  are  important 
here.  The  phase  velocity  is  defined  as 

oo 


which  in  this  case  is  1,  as  we  have  assumed  eo  —  po  —  c  —  1.  The  group  velocity  is  defined 
to  be 

C(k)  =  Vkw(k)  =  (7.5) 

It  is  a  well  known  fact,  that  the  propagation  of  energy  under  dispersive  partial  differential 
equations  is  governed  by  the  group  velocity  [35,  10,  36].  Asymptotically,  the  energy  asso¬ 
ciated  with  the  wave  vector  k  moves  at  the  group  speed  |C|,  which  in  the  present  case  is 
|  C  |  =  1.  Thus,  regardless  of  the  wave  number  k,  all  plane  waves  move  with  the  same  group 
speed  |  C  |  =  1. 


7.1  Dispersion  and  Stability  Analysis 

We  now  perform  a  similar  analysis  for  the  semi-discrete  system  (6.3),  where  we  consider 
exact  integration  in  time.  Let  us  assume  an  infinite  PML  in  the  region  x  >  0.  Thus,  <7^  =  0 
and  let  crx  =  a.  We  will  look  for  solutions  of  the  form 

V(x,y,t)  =  V(x,y)e'ut>  (7.6) 
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to  the  semi-discrete  system  (6.3).  Substituting  (7.6)  in  (6.3),  we  obtain  the  time  harmonic 
system 


i  uj(Bx,4>x 
iuj(Hx,ipx 


=  1C j 


1U 


1C V 


1  +  .  )  By,  Ipy  )  — 

liU  / 

\U(Hy,^y) 

i  +  -)  A 

1C JjJ 

ic o{E,<j>) 


dyE  ,  Ipa 

1  +  — 'j  Bx, 

iu/ 

9XE  ,  Tpy 
i  U(By,l/}y), 


(7.7) 


=  [Hx,dy 


Hy,  dx\ 


=  iw(D,</>). 

We  assume  that  a  is  a  piecewise  constant  function  of  x  with  jumps  at  x  —  Ih,  l  —  0, 1,  2, 
where  h  —  hx  —  hy  is  the  mesh  step  size.  Let 


<?i 


Value  of  a  on  ( Ih ,  (l  +  l)h),  if  l  >  0, 
0,  if  l<  0. 


(7.8) 


Using  the  definition  (3.3),  we  have 

sx  ,i  =  si  =  1  +  t — •  (7-9) 

1C 0 

Since  cry  =  0,  we  have  sy  =  0.  As  the  PML  is  in  the  half  space  x  >  0,  x  =  0  is  the  interface 
between  the  PML  and  the  interior  computation  region.  Let  us  define 


( 


\ 


MxUl,r 

MyUl^r 

SyUl^, 

Mzul>r 


4  ui  ,m  +  Ui- 1  + 

4  ui  +  Ui  ,771—  1  +  Ul,m+ 1, 

MyUt-  1/2  —  MyUl  +  l/  2  ,7715 

777+ 1/2  5 

MxMyUirri' 


(7.10) 


Consider  an  interior  super  element  as  shown  in  Figure  3.  Using  the  definitions  (7.10) 
in  (7.7),  we  obtain  the  following  system  of  equations  that  corresponds  to  the  space  discrete 
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Figure  3:  Dependency  diagram  for  an  interior  super  element. 


finite  element  scheme  (6.3): 


Mx  B{m+i  /  2 


MxHt  ,m+ 1/2 

SlMyBi+i/ 2 
MyHl+  1/2 
Sz  +  Sz-l 


r  -^x(-^Z,m+ 1  E/i  rn)) 

ujh 


Sl  +  Sz_l 


MxBinn-\- 1/2) 


—  i’m 

=  MyBl  +  l/  2 

,7715 


MzEitTn 


uh  v  y 
MzDir 


(7.11) 


Combining  the  equations  in  (7.11),  we  obtain  an  equation  in  if  by  eliminating  the  other 
variables,  which  is 

cu2/r2  /sz  +  sz_i 


6 


MzEim  — 


si  +  St-i 


(- MXE  l,m+l  2MxEim  +  MxEim_i) 


H - (MyEl+  1)TO  —  MyEl^jn) - (. MyEliln  —  MyEl- i,m).  (7.12) 

«Z  >Sz_i 
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Let  us  now  look  for  solutions  to  (7.12)  of  the  form 


Ehm  =  E^mh. 

Substituting  (7.13)  in  (7.12),  and  after  performing  some  algebra,  we  obtain 


OJh2 

6 


Sl\Sl  (4 Et  +  E{-i  +  Ei+i)  =  ~(Ei+ 1  -  Ei)  ~  —{Ei  ~  Et- 1) 

2  /  sz  Si-i 


where  the  coefficient  £  is  defined  as 


C  =  1 


12 


sm2(kyh/2) 


(7.13) 


(7.14) 


(7.15) 


uj2h2  \1  +  2  cos2(kyh/2) 

Let  kx  and  /c£ml  be  the  x  components  of  the  wave  vector  in  free  space  and  the  PML,  respec¬ 
tively.  Assuming  Ei  =  e~lkxhl  in  free  space,  substituting  in  (7.14),  with  a  =  0,  we  obtain  the 
dispersion  relation  in  free  space  to  be 

’ 2h 2  sin2(kxh/2)  |  sm2(kyh/2) 


UJ 


+ 


12  1  +  2  cos2(kxh/2)  1  +  2  cos2(kyh/2) 

The  dispersion  relation  for  the  FDTD  scheme  [33]  is 


(7.16) 


uj2h2 


=  siv?{kxh/2 )  +  siv?{kyh/2). 


(7.17) 


Thus,  depending  on  the  magnitude  and  direction  of  k,  the  numerically  computed  wave  has 
an  incorrect  phase.  As  a  result  a  plane  wave  of  the  form  (7.1)  will  generally  move  in  an 
incorrect  direction  at  an  incorrect  speed.  However,  from  (7.16)  we  can  observe  that  if  kxh 


and  kyh  are  small  we  have 


u+{kx,ky,h ) 


h 


2V3  (kxh)2  +  (kyhf 


12 


=  \  k2  +  k2  =  |k|,  as  h  — >  0 


(7.18) 


where  lv+  denotes  the  positive  solution.  Similarly  from  (7.17)  we  obtain  that  u+  — >• 
|k|,  as  h  0.  This  implies  that  the  effects  of  dispersion  can  be  reduced  to  any  desired 
level  if  we  choose  a  fine  enough  mesh.  To  derive  the  (angular)  frequency  00  for  the  fully 
discrete  scheme  (6.5),  we  observe  that  discretization  in  time  corresponds  to  replacing  uh  in 

(7.16)  by  2—  sin  ^  2  }  US  ^eno^e  ^ie  frecluency  f°r  bully  discrete  scheme  (6.5) 
by  a; At-  In  the  working  volume,  where  a  =  0,  we  get  the  dispersion  relation  to  be 

2 


h 


V3At 


sin 


/iuazAA 


sm2(kxh/2) 


+ 


sm2(kyh/2) 


1  +  2  cos2(kxh/2)  1  +  2  cos2(kyh/2) 


(7.19) 
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The  Courant  number  is  r)  =  cAt/h  (c  =  1).  Solving  for  u\t  (positive  solution)  in  (7.19)  we 
obtain 

,_i  fr]hoo+(kx ,  ky,  h) 


ul.  =  —  Sill 

Ai  At 


(7.20) 


where  u+  is  the  (positive)  solution  to  (7.16),  i.e., 


00  = 


2\/3  /  sin2  (kxh/ 2) 


+ 


sm2(kyh/2) 


h  V  1  +  2  cos2(kxh/2)  1  +  2  cos2(kyh/2) 


(7.21) 


We  note  that  rj  must  be  chosen  such  that  the  frequency  is  real.  This  in  turn  implies  that 
the  argument  of  sin'1  in  (7.20)  is  bounded  by  1.  We  thus  need 


V  (max(kxh,kyh)e[o,ir}x{o,ir]  \hu+(kx,ky,h)\ )  <  2 


2  1 

77  <  -==  =  -=  «  0.4082, 
v/24  V6 


(7.22) 


In  the  case  of  the  FDTD  scheme  a  similar  analysis  yields  the  stability  result 

1  cAt  1 

7?  <  —  ==►  —  <  —  «  0.7071, 
y/2  h  ~  ^2 


(7.23) 


with  c  —  1. 

Figure  4  compares  the  dispersion  in  the  FEM  scheme  and  the  FDTD  scheme  for  free 
space,  for  four  different  wave  propagation  angles,  9  —  0, 15,  30, 45  degrees,  with  77  =  0.4.  The 
x  axis  denotes  the  number  of  grid  points  per  wavelength,  L/h ,  where  L  is  the  free  space 
wavelength,  and  the  y  axis  is  the  numerical  phase  velocity,  vp,  normalized  to  the  speed  of 
light,  vp/c  (c  =  1).  We  note  that,  in  both  the  schemes,  the  phase  velocity  is  the  lowest  at  45 
degrees,  implying  that  the  dispersion  is  the  least  along  the  diagonal  of  the  mesh  elements; 
a  fact  that  will  be  evident  in  other  results  to  be  presented.  In  Figure  5,  similar  results  are 
presented  for  77  =  0.01.  We  note  that  the  dispersion  present  in  the  FEM  scheme  decreases 
as  the  value  of  77  is  decreased  from  0.4  to  0.01,  while  the  reverse  is  true  in  the  case  of  the 
FDTD  scheme.  The  major  effects  of  dispersion  are  seen  in  the  case  of  10  or  less  nodes  per 
wavelength.  In  all  cases,  as  L/h  becomes  large  the  convergence  of  vp  to  c  (=  1)  is  clearly 
seen. 

The  dispersion  relation  in  the  PML  is  obtained  in  a  similar  fashion.  Let  us  consider  0 


to  be  a  constant.  Assuming  E\  —  e  lk 
dispersion  relation  in  the  PML  to  be 


pml 


hi 


in  the  PML,  substituting  in  (7.14),  we  obtain  the 


ou 


2h2 


12 


sin2(£;Pml/i/2) 


+ 


sm2(kyh/2 ) 


s2  )  1  +  2cos2(&riV2)  l  +  2cos  2{kyh/2) 


(7.24) 
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Figure  4:  Comparison  of  the  dispersion  present  in  the  FEM  and  the  FDTD  scheme 
for  selected  angles  of  wave  propagation  with  rj  =  0.4. 
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Figure  5:  Comparison  of  the  dispersion  present  in  the  FEM  and  the  FDTD  scheme 
for  selected  angles  of  wave  propagation  with  r]  =  0.01. 


To  derive  the  dispersion  relation  for  the  fully  discrete  scheme,  we  again  replace  c oh  by 


h 

2—~  sm 
At 


f  uA  t 


to  get 


K2  .  2  ( uAtAt\  (  1  \ 
3(A tysm  {  2  J“Uv 


sin2(fcPmlh/2) 


+ 


sm2(kyh/2 ) 


1  +  2  cos2  ( kpxmlh/2)  1  +  2  cos2  (ky h/2) ' 


(7.25) 


7.2  Analysis  of  the  Phase  Error  and  Anisotropy 


The  phase  error  that  results  from  the  dispersion,  expressed  in  degrees  per  wavelength,  is 
defined  as 


360° 


k  —  k 
k 


(7.26) 


where  k  is  the  numerical  wave  number  with  which  the  plane  wave  propagates  in  the  numerical 
grid.  The  wave  number  for  the  continuous  case  is  k.  Figure  6  is  a  polar  graph  of  the  phase 
error  for  selected  values  of  L/h  =  2n/kh  as  a  function  of  9  for  the  FEM  and  FDTD  schemes, 
with  rj  =  0.4  (top)  and  q  =  0.01  (bottom).  From  Figure  6,  we  see  that,  for  both  the  schemes, 
the  smallest  error  occurs  when  the  plane  wave  traverses  the  elements  diagonally,  whereas 
the  largest  error  occurs  when  the  wave  is  propagating  along  an  axis  of  the  mesh  Similar 
observations  were  made  earlier  in  the  plots  of  vp/c  versus  L/h.  We  note  that  the  phase 
error  in  the  FEM  reduces  as  the  value  of  r /  is  decreased,  whereas,  in  the  FDTD  scheme,  the 
phase  error  increases  as  the  value  of  rj  is  decreased  from  0.4  to  0.01.  In  Figures  7  and  8,  this 
change  in  the  phase  error  for  the  two  schemes  is  more  evident.  In  Figure  7,  the  phase  error 
is  plotted  for  rj  =  0.4  (top)  and  rj  =  0.2  (bottom)  and  in  Figure  8  for  q  =  0.1  (top)  and 
rj  =  0.01  (bottom).  Again,  dispersion  effects  are  more  evident  for  the  case  of  10  nodes  per 
wavelength.  As  L/h  increases  the  phase  error  converges  to  zero,  for  all  angles  of  propagation. 
Thus,  the  effects  of  dispersion  can  be  reduced  to  any  desired  degree  by  considering  a  fine 
enough  mesh.  Similar  results  are  also  obtained  for  a  first  order  vector  finite  element  method 
applied  to  the  vector  Helmholtz  equation  [37]. 

Figure  9  is  a  log-log  graph  of  the  phase  error  5p,  versus  the  number  of  nodes  per  wavelength 
L/h,  for  selected  values  of  the  angle  6 ,  for  the  FDTD  and  the  FEM  schemes.  We  observe 
that,  for  both  schemes,  as  L/h  is  increased,  the  error  becomes  smaller.  For  large  values 
of  L/h  the  graphs  are  linear,  indicating  the  error  to  be  proportional  to  the  square  of  the 
inverse  of  the  number  of  nodes  per  wavelength,  i.e.,  to  {h/ L)2 .  This  implies  the  second  order 
convergence,  with  respect  to  (h/L),  of  k  to  k.  We  also  note  that  the  phase  error  is  lower  for 
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Figure  6:  Polar  graph  of  the  phase  error  in  degrees  per  wavelength  for  selected  values 
of  L/h,  with  r /  =  0.4  (top)  and  rj  =  0.01  (bottom). 

the  FDTD  scheme  than  for  the  FEM  scheme  when  rj  =  0.4,  but  as  rj  is  decreased  to  0.01, 
the  effects  of  dispersion  increase  for  the  FDTD  scheme. 

We  present  here  another  important  feature  of  numerical  dispersion,  which  is  the 
anisotropy  of  the  dispersion  with  respect  to  the  angle  of  incidence  of  the  propagating  wave. 
In  ordinary  wave  propagation,  energy  propagates  perpendicular  to  the  wave  front.  When 
there  is  anisotropy  dispersion,  the  angle  will  not  be  perpendicular.  The  anisotropy,  d,  is 
defined  as 

'0  =  max0<e<7r/25p  -  mino<0<^/2<5p-  (7.27) 

In  Tables  1  and  2,  we  present  the  maximum  and  minimum  values,  over  all  angles  of  prop¬ 
agation,  of  the  phase  error  5P,  for  10,  15,  20,  30,  and  40  nodes  per  wavelength.  We  also 
present  the  anisotropy  present  in  the  FEM  and  FDTD  schemes,  for  rj  =  0.4  and  rj  =  0.01, 
respectively.  In  Table  1,  with  rj  =  0.4,  we  see  that  the  maximum  and  the  minimum  values 
of  dp  are  larger  in  the  case  of  the  FEM  scheme;  however,  the  anisotropy  is  less  in  the  case  of 
the  FEM. 

In  Table  2,  with  rj  =  0.01,  we  see  that  the  maximum  and  minimum  values  of  the  FDTD 
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Figure  7:  Phase  error  in  degrees  per  wavelength  as  a  function  of  the  angle  6  for 
selected  values  of  L/h ,  with  r/  =  0.4  (left  box)  and  rj  =  0.2  (right  box). 


L/h 

FEM 

FDTD 

Max 

Min 

d 

Max 

Min 

d 

10 

6.5733 

3.8228 

2.7505 

5.2043 

2.0599 

3.1444 

15 

2.9930 

1.7204 

1.2726 

2.2548 

0.9044 

1.3503 

20 

1.6981 

0.9720 

0.7261 

1.2573 

0.5066 

0.7507 

30 

0.7594 

0.4334 

0.3260 

0.5539 

0.2245 

0.3309 

40 

0.4281 

0.2441 

0.1841 

0.3117 

0.1261 

0.1856 

Table  1:  Anisotropy  for  rj  =  0.4,  for  selected  values  of  L/h. 


scheme  are  now  larger  than  the  respective  values  in  the  FEM  scheme.  The  anisotropy  is  also 
smaller  for  the  case  of  the  FEM. 

7.3  Calculation  of  the  Reflection  Coefficient 

In  this  section  we  study  the  properties  of  the  discrete  FEM-PML  model  by  performing  a  plane 
wave  analysis  to  calculate  the  reflection  coefficient.  In  the  discrete  setting  the  PML  model  is 
no  longer  perfectly  matched  since  the  discretization  introduces  some  error  which  manifests 
itself  as  spurious  reflections.  There  is  also  error  that  is  introduced  due  to  the  termination  of 
the  PML.  We  study  the  errors  introduced  in  the  discrete  model  by  calculating  the  reflection 
coefficient  of  an  infinite  PML  (to  study  the  errors  caused  by  the  discretization)  as  well  as 
the  reflection  coefficient  of  a  finite  PML  (to  study  the  errors  introduced  by  terminating  the 
PML). 
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Figure  8:  Phase  error  in  degrees  per  wavelength  as  a  function  of  the  angle  6  for 
selected  values  of  L/h ,  with  77  =  0.1  (left  box)  and  77  =  0.01  (right  box). 


L/h 

FEM 

FDTD 

Max 

Min 

7? 

Max 

Min 

7? 

10 

5.6705 

2.8991 

2.7714 

6.2007 

3.0302 

3.1705 

15 

2.5812 

1.3043 

1.2769 

2.6850 

1.3298 

1.3552 

20 

1.4643 

0.7368 

0.7275 

1.4970 

0.7447 

0.7522 

30 

0.6548 

0.3285 

0.3263 

0.6612 

0.3299 

0.3312 

40 

0.3691 

0.1850 

0.1841 

0.3711 

0.1854 

0.1857 

Table  2:  Anisotropy  for  77  =  0.01,  for  selected  values  of  L/h. 


For  simplicity  we  again  assume  in  this  section  that  eo  =  Ho  —  1.  Let  us  also  assume  an 
infinite  PML  in  the  region  x  >  0.  Thus,  ay  =  0  and  let  crx  =  er.  To  calculate  the  reflection 
coefficient  for  the  infinite  PML,  we  look  for  solutions  to  (7.14)  of  the  form 


e~lkxhl  +  Reik*M,  for  l  <  0, 
Te-fe"mlw,  for  l  >  0. 


(7.28) 


The  reflection  coefficient  is  R,  and  T  is  the  transmission  coefficient.  Consider  the  equations 
associated  to  the  node  at  the  interface  l  =  0  and  one  node  each  on  either  side  of  the  interface 
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Figure  9:  Log-Log  plot  of  the  phase  error  in  degrees  per  wavelength  as  a  function 
of  L/h  for  selected  values  of  the  angle  of  incidence  6  with  r)  =  0.4  (left  box)  and 
rj  =  0.01  (right  box). 


at  l  —  1,  and  l  =  —  1.  From  (7.14)  we  have 

_^I(4jF_1  +  E—2  +  E0)  =  (4  -  £U)  -  (E-i  -  E_2), 

o 

_(urfr  /«o\  ^  +  ^  +  ^  =  i(4  _  ^  ^  _  44  (7.29) 

o  V  2  /  So 

_CcWr  /  s_i+^o\  (44  +  4  +  4)  =  1(4  -  4)  -  1(4  -  Eo), 

6  V  2  /  «i  s0 

where  £  is  defined  in  (7.15),  and  si  is  dehned  in  (7.9).  Substituting  for  4  from  (7.28)  in 
(7.29)  we  obtain  three  equations  in  the  unknowns  4>  R  and  T.  Solving  these  resulting 
equations  for  R ,  we  can  show  that  the  reflection  coefficient  has  the  Taylor  series  expansion 

R  =  — — l^-(co»2  -  kl)a(a  +  2iuj)h2  +  - -*—0-2(<t  +  2icu)(u;2  -  k2)3/2h3  +  0(hA).  (7.30) 

loc dz  y  48cL  y 

The  formula  (7.30)  implies  that  the  reflection  coefficient  is  proportional  to  h2.  Thus,  the 
discrete  PML  model  is  consistent  with  the  continuous  model. 

Next,  we  study  the  effects  of  terminating  the  PML  by  a  PEC.  This  amounts  to  setting 
E  —  0  at  the  boundary  x  =  5  =  Mh  of  the  PML,  i.e.  Em  —  0.  To  obtain  the  reflection 
coefficient  we  write  equation  (7.14)  for  all  the  nodes  in  the  PML  as  well  as  for  the  node  at 
the  interface  of  the  working  volume  and  PML,  4,  and  node  E-\  in  the  working  volume 
which  is  h  distance  away  from  the  interface.  Assuming  that  we  know  the  value  of  EL 2  we 
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obtain  a  system  of  equations 


AE  =  ~(u2h2(  +  6)E_2e1. 


(7.31) 


In  the  above  E  =  [E-i,  E(h  E±,  •  •  •  ,  Em_i]t,  eq  =  [1,  0,  0,  •  •  •  ]T  and  the  matrix  of  coefficients 
obtained  from  (7.14)  is 


where 


A  = 


6-i  c0  0 

a_  i  60  ci  0 

0  a0  b\  c2 


0  am- 2 


a,  =  |  (  ftti-tq  +  1 

2  sz 


6Z=  [4a;2h2'  ' Sl"1  +  Si 


-6  -  + 


Sl  <sz_i 


I  2r2  /  Sl-1  +  Sl—2  6 

Q  =  w  h  |  - x - h 


(7.32) 


(7.33) 


2  si-ii 

We  can  solve  system  (7.32)  for  the  value  of  R  by  using  (7.28)  for  l  =  —1  and  l  =  —2.  In  this 
case  the  absolute  value  of  the  reflection  coefficient  is  calculated  to  be 

1  +  (u2h2(  +  6)k  eik*h 


\R\  = 


1  +  (u)2h2(  +  6  )/ce  1 


— i  kTh 


(7.34) 


where  k  is  the  hrst  diagonal  entry  in  A -1. 

Figure  10  plots  the  reflection  coefficient  in  decibels,  Db  (i.e.,  20  log10  /?),  versus  the 
number  of  nodes  per  wavelength  L/h  for  different  values  of  Rq,  the  reflection  coefficient 
at  normal  incidence.  Figure  11  plots  the  reflection  coefficient  in  Db  versus  the  angle  of 
incidence  6.  In  these  figures  we  compare  the  reflection  coefficient  for  the  FEM  scheme  with 
the  reflection  coefficient  for  the  TE  mode  of  the  Zhao-Cangellaris’s  PML  model  using  the 
FDTD  scheme  which  was  presented  in  [14]. 

We  note  that  the  numerical  reflection  coefficient  converges  to  the  reflection  coefficient  of 
the  continuous  model,  which  is  Rq0s9  as  we  increase  the  value  of  N.  We  also  note  that  as  6 
approaches  the  value  tt  the  numerical  reflection  coefficient  approaches  the  value  1.  This  is 
a  well  known  behavior  of  PML  models,  i.e.,  waves  that  are  propagating  transversely  to  the 
interface  between  the  domain  of  interest  and  a  single  PML,  are  not  absorbed  by  the  PML. 
However,  these  waves  get  absorbed  into  the  corner  regions  where  two  PML’s  overlap.  The 
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Reflection  Coefficient  R  in  Db 


Numerical  Reflection  Coefficient  for  0  =  0  Numerical  Reflection  Coefficient  for  0  =  jc/4 


Number  of  points  per  wavelength,  N  Number  of  points  per  wavelength,  N 

Figure  10:  Numerical  reflection  coefficient  for  6  =  0  (left)  and  9  —  7t/4  (right). 
We  note  that,  as  we  increase  the  number  of  nodes  per  wavelength,  the  numerical 
reflection  coefficient  in  both  cases  approaches  Rq. 


Numerical  Reflection  Coefficient  for  R0  =  1  .Oe-2  Numerical  Reflection  Coefficient  for  R0  =  1 .0e-4 


Angle  of  Incidence  in  Radians  Angle  of  Incidence  in  Radians 


Figure  11:  Numerical  reflection  coefficient  for  i?o  =  10  2  (left),  Rq  =  10  4  (right). 
As  L/h  is  increased,  the  numerical  reflection  coefficient  converges  to  Rq°s6  (Exact). 
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plots  were  obtained  by  considering  PML’s  that  are  one  wavelength  thick,  i.e.,  the  number 
of  nodes  per  wavelength  is  the  number  of  nodes  in  the  PML.  The  polynomial  grading  (3.4) 
was  chosen  for  o  with  rn  —  2  and  <rmax  as  in  (3.5)  with  Z  —  1.  The  PML  is  in  the  region 
x  >  0.  Thus,  x0  =  0  in  (3.4)  for  a  =  x. 

The  power  of  the  polynomial  grading  for  a  has  an  effect  on  the  numerical  reflection  errors 
as  is  known  from  analysis  of  PML  models  discretized  using  the  FDTD  scheme.  In  Figure  12 
(top),  we  plot  the  reflection  coefficient  for  the  FDTD  and  FEM  schemes  against  the  number 
of  points  per  wavelength  for  m  =  2.55  and  m  =  8.65,  with  6  =  0  degrees,  Rq  =  1CT4 
in  each  case.  In  Figure  12  (bottom),  we  plot  the  reflection  coefficient  (FEM)  against  the 
number  of  points  per  wavelength  and  m  for  6  =  0  degrees,  Rq  =  10-4.  Both  these  figures 
demonstrate  the  strong  influence  of  the  power  m  of  the  polynomial  grading  for  a  on  the 
reflection  coefficient. 


8  Absorption  of  a  Pulse  on  the  Boundaries  of  a  Com¬ 
putational  Domain 

The  numerical  experiment  described  in  this  section  evaluates  the  performance  of  the  UPML 
when  a  pulse  strikes  the  boundaries  of  a  computational  domain.  We  measure  the  amount 
of  reflection  that  an  outward  propagating  pulse  produces  as  it  moves  from  free  space  to  a 
boundary  surrounded  by  absorbing  PML’s  as  in  Figure  1. 

We  choose  our  domain  D  to  be  the  square  [0,12]  x  [0,12],  with  a  source  located  at 
the  center  (6,  6)  of  the  square.  The  domain  is  surrounded  by  absorbing  layers  on  all  four 
boundaries.  We  discretize  the  problem  with  a  rectangular  grid  composed  of  90  x  90  square 
elements  of  step  size  h  =  2/15  and  the  time  step  is  At  =  0.04/c,  which  satisfies  the  CFL 
condition  (7.22).  The  source  is  taken  to  be  the  function  [12], 


f(x,y,t )  =  fi(x,y)f2(t), 


where 


In  the  above,  /0 
defined  as 


h(t)  = 


-2Tr2f^(t-tQ)e-”a^t-^\  iff  <2f0, 
0,  if  t  >  2 10- 


(8.1) 


20  h 


is  the  central  frequency  and  t0  =  1  / f0.  The  function  fi(x,y)  is 


h  (x,y)  =  e-VP-6)2+D-6)2_ 


(8.2) 
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Numerical  Reflection  Coefficient  for  Normal  Incidence 


Number  of  points  per  wavelength,  N 


Figure  12:  (top)  Numerical  reflection  coefficient  for  FDTD  and  FEM  in  two  cases; 
m  =  2.55  and  m  =  8.65.  (bottom)  Numerical  reflection  coefficient  for  different 
points  per  wavelength  and  different  m  values.  In  both  figures  9  —  0,  Ro  —  10-4. 
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Figure  13:  Comparison  of  the  L2  error  for  the  UPML  with  a  mixed  finite  element 
scheme  and  the  split  field  PML  with  the  FDTD  scheme  on  a  90  x  90  cells  grid  (left) 
and  for  a  180  x  180  cells  grid  (right).  We  note  that  as  the  grid  is  refined  and/or  as 
the  number  of  PML  cells  is  increased  the  error  in  the  FEM  and  FDTD  schemes  is 
identical. 

We  obtain  a  reference  solution  by  using  a  similar  finite  element  scheme  for  the  TM  mode 
of  Maxwell’s  equations  on  a  larger  domain  f 1r  containing  360  x  360  square  elements,  and 
the  same  mesh  step  size  and  time  step.  The  domain  Dr  is  terminated  using  PEC  conditions 
on  its  boundary.  We  have  used  the  polynomial  grading  (3.4)  for  o  with  the  optimal  value  of 
cr max  as  given  in  (3.6)  with  m  =  3.5. 

The  L2  norm  of  the  error  clue  to  numerical  reflections,  which  arise  clue  to  the  finite  PML 
terminated  by  PEC  conditions,  is  obtained  by  subtracting  at  each  time  step  the  field  E  at 
any  grid  point  inside  D,  from  the  field  E  at  the  corresponding  point  in  Dr,  taking  the  square 
of  this  difference  and  summing  such  differences  over  all  grid  points  in  D.  We  do  the  above 
for  three  PML’s  containing  4,  8  and  16  cells.  A  comparison  is  presented  with  respect  to  the 
split  field  PML  (SF)  of  Berenger,  using  the  same  test  problem.  The  reference  solution  for 
the  split-field  case  is  constructed  in  a  similar  way.  Figure  13  shows  the  L2  error  between  the 
two  reference  solutions  (Reference  Error)  for  the  split-field  and  the  finite  element  schemes, 
and  the  L2  error  of  the  two  schemes  for  4,  8  and  16  cell  PML’s  each. 

From  Figure  13  (left),  we  can  see  that  the  reference  error  (discretization  error)  dominates 
for  about  250  time  steps.  After  this,  as  the  wave  exits  the  computational  domain,  the 
reflection  error  due  to  the  PEC  backed  PML  takes  over.  We  have  used  20  nodes/ wavelength 
(i.e.,  L/h  =  20)  in  our  calculations.  As  can  be  seen  for  a  16  cell  PML  the  reflection  error 
is  lower  than  the  reference  error.  Figure  13  (right)  shows  the  L2  error  between  the  two 
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Figure  14:  Propagation  of  the  wave  front.  A  top  view  of  the  solution  is  plotted  at 
different  time  steps.  The  magnitude  of  the  solution  is  the  same  for  all  time  steps. 


reference  solutions  (Reference  error)  of  the  split-field  and  the  finite  element  scheme,  and  the 
L2  error  of  the  two  schemes  for  4,  8  and  16  cell  PML’s,  for  a  refined  discretization.  In  this 
case,  h  =  1/15  and  At  =  0.02/c.  From  Figure  13  (right)  we  can  see  that  a  four  cell  PML 
provides  a  good  absorbing  layer. 

In  Figure  14,  we  plot  the  propagation  of  a  pulse  on  a  180  x  180  cells  domain  backed  by 
an  eight  cell  PML  obtained  using  the  FEM  scheme.  The  wave  front  completely  disappears 
from  the  domain,  as  seen  in  the  subplot  corresponding  to  320  time  steps.  All  subplots  are 
plotted  at  the  same  magnitude. 

9  Conclusion 

In  this  paper,  we  have  presented  and  analyzed  a  mixed  finite  element  scheme  for  the  numerical 
solution  of  the  2D  TM  mode  of  the  uniaxial  PML.  We  have  stated  energy  estimates  under 
certain  assumptions.  Based  on  these  estimates  and  the  dispersion  and  reflection  coefficient 
analysis  of  the  scheme  we  can  conclude  that  the  proposed  scheme  has  absorbing  properties 
that  are  comparable  to  those  of  PML  models  discretized  using  FDTD.  The  extension  of 
the  mixed  FEM  to  3D  is  straightforward  and  uses  a  combination  of  Nedelec’s  elements  and 
Nedelec-Raviart-Thomas  elements  for  the  discretization  of  the  electric  and  magnetic  fields, 
respectively. 
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As  in  the  case  of  the  FDTD  method,  we  also  find  here  that  the  choice  of  the  PML 
conductivity  affects  the  numerical  reflection  errors.  Thus,  a  more  rigorous  analysis  is  needed 
to  determine  the  optimal  choice  of  the  polynomial  approximation  to  the  PML  conductivity 
that  will  minimize  the  numerical  reflection  errors  that  are  generated. 
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